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/3. The latter secular coefficients are related to the traces of powers of U by Newton's well- 
known formulae. While the traces tend to have Gaussian distributions and to be statistically 
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coefficients exhibit strong mutual correlations due to Newton's mixing of traces to coefficients. 
These results might become relevant for current efforts at combining semiclassics and random- 
matrix theory in quantum treatments of classically chaotic dynamics. 
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1. Introduction 

Circular ensembles of unitary matrices were first considered by Dyson B and described in detail 
by Mehta |^] . They are used to describe the quantum statistics of periodically driven systems (see 
H and references therein) and of scattering processes [Q . It is the unitary Floquet operator in the 
first case and the unitary S'-matrix in the second that the random unitary matrices in question 
attempt to mimic with respect to certain more or less universal properties. Systems with global 
chaos in their classical limit and with time reversal invariance either present or strongly broken 
exhibit the greatest degree of universality in their statistical properties and tend to fall in one of 
the universality classes represented by the circular orthogonal, symplectic, or unitary ensembles 
(COE, CSE, CUE). The so-called Poissonian ensemble (CPE) of diagonal unitary matrices with 
independent unimodular eigenvalues has also found applications for certain classically integrable 
systems. There is a recent interest in the analysis of intermediate ensembles of unitary matrices 
which describe cross-over between different universality classes || , |J , . Direct links between 
corresponding Gaussian and circular ensembles have been established lately j^], j|. 

Unitarity constraints imposed on a random matrix of small size N cause a significantly 
non-Gaussian character of the distribution of matrix elements (To), Jl^]. Moreover, various 
statistics of unitary matrices from the different circular ensembles depend strongly on the matrix 
size, and for small N differ a great deal from the asymptotic large N properties. This is 
in contrast to the Gaussian ensembles of Hermitian matrices j|] which tend to display lesser 
sensitivity to the matrix dimension. 

For example, the normalized level spacing distribution P(S) suffers a cut-off at S = N, and 
is close to the Wigner-like distributions only for TV upwards of roughly 10. Moreover, Baranger 
and Mello have shown [jl3| that the distribution of transmission intensities P{T) is non-Gaussian 
for N = 2 and 4, while coming close to Gaussian for larger N. 

We were led to study the secular polynomials of random unitary matrices in an attempt at 
constructing semiclassical quasienergy spectra for the kicked top under conditions of classical 
chaos. Periodic-orbit theory can be invoked only insofar as periodic orbits are available. In 
practice, like for any system for which no simple symbolic dynamics is known, one can hope 
to find all periodic orbits with periods no longer than maybe n m ax ~ 10, due to the infamous 
exponential proliferation. These would allow to semiclassically evaluate traces of powers of the 
Floquet operator, trF n , with n up to n m ax- A Hilbert space of dimension N = 2n m ax is 
then accessible since Newton's formulae [H (see Section |f) allow to express the first Umax 
coefficients of the secular polynomial in terms of the first n m ax traces; the so-called self- 
inversiveness of the secular polynomials of unitary matrices EE) then yields the second half 
of the set of coefficients. The practical applicability of periodic-orbit theory would thus seem 
severely limited: Inasmuch as the dimension N of the Hilbert space is proportional to the 
effective size of quantum fluctuations (formally, N oc 1/h) one runs out of periodic orbits just 
when the semiclassical approximation begins to have a raison d'etre. 

To ease the dilemma just described one might hope to increase the size of the Hilbert space 
by throwing dice, according to random-matrix theory, for traces with exponents n m ax < n < 
N — rimax- To prepare for such a "marriage" of semiclassical approximations with random- 
matrix theory we here propose to study secular polynomials of random unitary matrices from 
various ensembles. A previous first step in this direction was taken in ]l6| , where the means 
and the mean squares of the coefficients mentioned were calculated for the CPE and the CUE. 
A related study was presented by Bogomolny, Bohigas, and Leboeuf fll5|] , fl7f who gave the 
distribution of the roots of random self-inverse polynomials. 

Even though the bulk of the work to be presented is analytical we performed extensive 
comparisons with numerical data on ensembles of unitary matrices. To generate the data we 
constructed COE and CUE matrices with the algorithm given in |1| and CSE matrices as 
described in |l9| . 



2. Coefficients of secular polynomials 
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2.1. Theory 

The secular polynomial of a unitary matrix U of size N is defined as 



N 



det(U — A) = ]T(-A)" aA r_„ = H (e^ - A) (1) 

n— i 

where the are the eigenphases of U and a = 1. We are interested in the statistics of the 
secular coefficients a n due to various ensembles of random matrices. We shall characterize these 
ensembles by their joint densities of the N eigenphases tpi, 

dP N {(p u <p 2 ,... <p N ) = C(J3, N) 11 |e^ - e'w f, (2) 

i<j 

where C(f3,N) is a normalization constant while /3, the so-called degree of level repulsion, 
distinguishes the ensembles: j3 = for the Poissonian ensemble for which the eigenphases 
are independently and uniformly distributed over the interval [0,27r); the circular orthogonal, 
unitary, and symplectic ensembles are characterized by, respectively, (3 = 1, /3 = 2, and (3 = 4. 
Occasionally we shall allow the parameter (3 to range freely among the real numbers. All of the 
ensembles in consideration are homogeneous inasmuch they do not distinguish any particular 
value of any eigenphase. It follows immediately that the ensemble means of all coefficients of 
the secular polynomial vanish, 

dcpi . . . d(p N a,Tfdfr(<pi,ip2, . . . <pn) = 0. (3) 



Next, we are interested in the ensemble averages a n a* m which can be obtained from the 
generating function 



= det(U -\){U- /i)t = Y[(c iVk - A)(c i v fc - fx)* (4) 

k 

= f d<p 1 ...d<p N dfa<pi,<p 3) ...<p N )W&>-\)(e i *>-tf. 

For A = e lv ,/j, = e lx it follows from the periodicity of the integrand that the integral depends 
on the phases tp, x only through the variable x = e 1 ^ _x - ) . Thus our generating function may be 
written as 

JV 

pn&) = n^ 1 - - e_ix ) = n ffa> x ) = e xn \^ ( 5 ) 

i i n— 

where we have introduced the auxiliary function 

f(^x) = (e^-x)(e-^-l) . (6) 



Hence all correlations a„a^ for m ^ n vanish. Moreover, one easily shows 



Kl 2 = \a N -n\ 2 (7) 

which is in accord with the so-called self-inversiveness, ajv_„ = a^apf [[l4j; being a consequence 
of but slightly weaker than unitarity self-inversiveness entails each root of a polynomial to either 
lie on the unit circle of the complex plane or to be accompanied by its inverse as another root. 
The variances \a n \ 2 are most easily calculated in the Poissonian case (3 = 0, 

,2tt n N 



P° N (x) = C(0,N) 



d<pf((fi,x) 



(l+xf. (8) 



Therefore \a n \ 2 = (^) for f3 = (T§. 

Next, we turn to the case (3 = 2 which was already treated in reference [Q. Observing that 
the function FL/CVii 3 ') which we want to average is symmetric in the N phases tpi and that 
the density (P N may be written as a product of two Vandermonde determinants we have Q , S 

pUx) = L jkw det (e i(m - iv - n ( 9 ) 
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with m,n = 1,2, . . .N labelling the rows and columns of the determinant. The integral over 
ip m can now be pulled into the mth row of that determinant whereupon we immediately get 

N 

Pjj{x) = det ((1 + x)S(m -n)- S(m -n+1)- x5(m - n - 1)) = ^ x N (10) 

where 5(m — n) denotes the Kronecker delta. The foregoing generating function entails the 
variance \a n \ 2 = 1 for p = 2. 

A little more effort is required for the orthogonal and symplectic cases (/3 = 1,4). In the 
orthogonal case we again employ the symmetry of the integrand of the TV-fold integral in (||) to 
rewrite that integral, which goes over the hypercube < tpx, if2, ■ ■ • , <Pn < 27r, as ./V! times one 
over the hypertriangle 2n > <pi > p>2 > ■ ■ ■ > <pn > 0. Within the hypertriangle the product of 
differences takes the form 

Y[ \ e iVk -e iv '\ = ]]_2sm f^y^) = i i ^ i ^ det (e imifil , . . . ,e imVN ) ; (11) 

k<i k<i ^ ' 

the second member in this chain of equated expressions is manifestly positive whereupon the 
modulus operation can be dropped; the determinant in the last member of (|ll|) has its rows 
labeled by the parameter m which runs in integer steps between and • the label for the 
columns is the one on the N integration variables. We may thus write the generating function 
as 

Pk(x)=NlC(l,N) J d^ 1 ...d^i^ i2 det(e i "^,...,e im ^)[]/(^,x). (12) 

tpi><P2—><PN fc 

For the symplectic case we find it expedient to extend the TV-fold integral in (||) to a {2N)~ 
fold one, 

1...2JV / / 

/ ip k - (pi 



p4( x ) = c(4,7V)iV! [ d^L.-dyjajv f[ U sin 

k<* P ^ 

V1>V2>--->¥'2JV 



V 2 



(13) 



N 

■ 

k=l 



n 



f(f2k,x) ( q — — ) S {(f2k-l ~ V2fc ~ e) 



2k— 1 



where e is to be sent towards zero from above. To see the equivalence of the foregoing expression 
to Pff as given by (||) with P = 4 we simply integrate by parts with respect to the phases with 
odd labels. Only those terms survive for which the TV differentiations have turned precisely those 
TV sine functions into cosines which are assigned vanishing arguments by the delta functions; 
the remaining sine functions then come in quadruples like, symbolically, S13S14S23S24 — * s| 4 . 
We should note in passing that the delta functions in the foregoing representation of P^ reflect 
Kramer's degeneracy and that in our definition of the secular polynomial for the symplectic 
ensemble each of the 2-fold degenerate eigenphases appears only once. 

The integral representations ( [l2| ) and (13) are convenient starting points for an explicit 
evaluation of the generating functions Pj^ and P^ . We propose to start with the slightly easier 
symplectic case. 

In order to actually evaluate the 2n-fold integral in (13) we once more employ the identity 
( |TT| ) as extended to 2iV ordered phases; we then integrate over every second phase to get rid 
of the delta functions. Upon renaming the remaining integration variables as (f2k ~^ fk and 
letting the positive infinitesimal e go to zero we arrive at 

Ph(x) = (-1) N C(4,N) / d N V det(me im ^\e m ^\me im ^,e im ^,..)T\f{ip k ,x) (14) 
Jo fe=i 

with |m| < (2N — l)/2. Note that we have returned to an integral over an N dimensional 
hypercube, exploiting the symmetry of the integrand. The difficulty to be coped with now lies 
in the fact that each of the N integration variables appears in a pair of rows of the 2N x 2N 
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determinant. At this point it is helpful to express that determinant by a Gaussian integral over 
complex Grassmann variables as 

det M = f (jl d V*k^ cxp (- v* k M klV )j . (15) 
Accounting for the matrix M from ( |l4| ) and doing the 2JV integrals over the rf we arrive at 



d^VdetAf = 




exp 



is 

m.m 1 



(16) 



with the antisymmetric 2JV x 2JV matrix 

271 

A mm , = {m! -m) J &ip^ m+m ' )v f{v,x) 


= 27r(m' — to) ((1 + x)S(m + m ) — <5(m — to' + 1) — xS(m + m — 1)) 



(17) 



the labelling of the rows and columns of A is inherited from the ordered-phases form of the 
Vandermonde determinant ([ll]): Both to and ml run in integer steps from —(27V — l)/2 to 
(27V — l)/2. The remaining Gaussian integral in (16) is easily recognized as the Pfaffian v 7 det A 
of the antisymmetric matrix A, such that the generating function in search takes the form 



P*(x) = C(4, JV)JV! VdetA. 



(18) 

The sign of the Pfaffian Vdet A must be chosen such that Pj^(x) is positive for positive values 
of x, according to the definition of the generating function. 

In a first attempt at evaluating the Pfaffian one may rejoice in the ease in finding it for small 
values of N which latter suggest a surmise for the general variance, 

'N\ 1 • 3 . . . (2n - 1) 

;/ 1 (_\\ 1 wis - H)...(2N -2n + l)' 



12 - 



for/3 



(19) 



(2N- 1)(2AT-3)...(2AT- 

This conjecture will be proven in the appendix. 

We finally turn to the orthogonal case, taking up the integral representation jl^ ) for P^(x). 
For the sake of simplicity let us assume an even dimension N. We start with integrating over 
the angles ipi with even indices i, pulling the (2A;)th such integral into the (2fc)th row of the 
determinant; while that integral at first appears as going over the interval c/?2fc+i < ^2k < <^2fc-i 
we can hurry to replace the lower limit with 0, simply by adding the iVth column to the (N— 2)th, 
the resulting (N — 2)th to the (JV — 4)th and so forth and thus obtain 



P^(a;) = C(l,JV)JV! 



d(^id(^3 . . . dtpN-i i 2 



(20) 



ipi>ip3>...>tfN-l 

■ <!«'! i . ' i, ^\/(^i,^),|d^e im ^/(^ a ;),e im ^/^3,^,y'd^e i ^/(^,x),. 

Now the integrand is symmetric in the remaining N integration variables whereupon we may 
extend the integration range to the N dimensional hypercube of edge length 2n and make up 
by dropping the factor JV!, as we had previously done in ([lj); the analogy with ( |l4| ) in fact goes 
much further: Once more, every integration variable appears in two rows of the determinant. 
Going through precisely the same reasoning as before we again incur a Pfaffian form, 

(21) 



P^(x) = C(l,N)N\ VdetB 



with the slightly more unpleasant antisymmetric matrix 



B,, 



(22) 



() 



here, the sign function ensures the antisymmetry of the matrix B with \m\ < (N — l)/2. Again, 
the Pfaffian suggests, by its easily evaluated form for small dimensions N, a guess for the 
variances, 

^- F=i+ n^_n) for/3 = 1 (23) 

We refer to the appendix for the proof of that surmise for N even or odd. 

Upon inspecting the variances found above for (3 = 0, 2 and conjectured for (3 = 1, 4 we were 
led to extrapolate to arbitrary non-negative (3 as 



12 - 



1 • (1 + f)(l + /?)... (l + i^) 



(l + (N l)f) (l + (N - 2)§ ) . . . (l + (N - n)f ) 



N \ r(n+|)r(jV-n + |) 

rj r(|)r(jv + |) 



(24) 



This expression has poles at /3 = , > • • • aT^TT an( ^ zeros at/3=^p,^,^,... — j for 
ft < It is thus analytic and positive for all positive (3, and it goes to zero for (3 — > oo. This 
looks like strong evidence for the general validity claimed before. We shall actually turn the 
conjecture into a theorem in the appendix. The proof will be based on the fact that (24) is 
equivalent to the differential equation 



for the generating function P^(x) 



2.2. Numerical results 

We constructed random unitary matrices U of different sizes according to the algorithm 
developed in |Q for the CUE and the COE and later generalized for the CSE JTsQ . For each 
such matrix we calculated a complete set of N secular coefficients a n by first computing the 
traces of arbitrary powers, t n — tr(f/ n ), via either matrix multiplication or diagonalization; 
Newton's formulae O] (see next section) then led to the a n . 

Precise estimates of the variance of any random variable require a much larger sample than 
estimates of the mean. We therefore present numerical results obtained for large samples of 
relatively small random matrices (N ~ 20), although some computations done for N ~ 200 
provide similar results. 
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N 

Figure 1. Variance of the distribution of the first coefficient |ai| 2 as a function of the matrix 
size N for circular ensembles: COE (□), CUE (O ) and CSE (A). Solid lines represent analytical 
results and dotted lines the asymptotic behaviour. 

Our above formulae for the variance \a n \ 2 of the secular coefficients involve the index n of 
the coefficient and the matrix size N. Figure [l] represents the variance of the first coefficient as 
a function of N, while Figure || shows the dependence on n for a fixed matrix size. Due to the 
property (0) of self-inversiveness the latter curve is symmetric about n = N/2. 
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Figure 2. Variance of the distribution of the nth coefficient |a n | 2 obi™!™ fr< 
of size TV = 20 typical of COE (□), CUE (o ) and CSE (A). 



The data of all three ensembles coincide (up to a statistical errosi aaith 
predictions. Note that the width of the distribution of coefficients deilielfes vith increasing 
degree of repulsion and is smallest for the symplectic ensemble. Additi mal numerical 
investigations confirmed the expectation that the phases of the a n are distrit uted 
in the range [0, 2ir) for any canonical ensemble. 



o 



b) 



"m 40000 matrices 



the theoretical 




Figure 3. Variance of coefficients for cross-over CPE - CUE: jai| 2 (a) and | q.2 | 2 (b) as functions 
of the matrix size N. Control parameter <5 is equal to 0.0 (Poisson, O), 0.2 (#), 0.4, 0.7 I and 
1.0 (CUE, O ), while the level-repulsion parameter j3 fitted simultaneously for both coefficients 
equals 0.0, 0.06, 0.24, 0.69, and 2.0, respectively. 



Even though the formula (24) for the variance is primarily meant to cover the four canonical 
ensembles /3 = (CPE), 1 (COE), 2 (CUE), and 4 (CSE), we could not resist the temptation 
to test its usefulness for intermediate cases. To this end we constructed an ensemble of 
unitary matrices interpolating between the Poisson and unitary ensembles according the method 
presented in [^p| . This intermediate ensemble depends on one control parameter 5, varying from 
(CPE) to 1 (CUE). Figure | shows the dependence of the variance of the first two coefficients on 
the matrix size N for S equal to 0.0, 0.2, 0.4, 0.7, and 1.0. For each case the value of the parameter 
(3 chosen to fit the N dependence of |ai| 2 , inserted into (24) provides a fair approximation for 
|a.2 1 2 , and similarly for subsequent coefficients. This astonishing fact reveals a certain validity 
of the general formula (24) with non-integer values of (3 for ensembles in between the usual 
universality classes. 



3. Traces of powers of matrices from circular ensembles 

The characteristic polynomial of a matrix U is related to the traces of its powers 



tr(ET 



N 

) = £ einV * 
i=l 



(26) 
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by 



det(U — A) = (-A) iv exp ( trln(l - -£/) 



(27) 



Expanding both sides in powers of A one finds the explicit relations between the a n and the t n 
which were already established by Newton mM. A compact representation is 



(28) 





h 


i 








■ 


1 

a n = —: 
n\ 




ti 


2 





• 


h 


t 2 


h 


3 


■ 








tn-2 


tn-3 


• h 



We infer from (|27| ) that there are only AT independent traces t±, . . . ,tisr. Moreover, for unitary 
matrices for which the eigenphases <pi are real the number of independent complex parameters 
is again reduced by a factor 1/2: The first N/2 traces suffice to determine all N coefficients a n 
(cf©). 

Clearly, the t n all vanish in the mean for all of the circular ensembles considered here, due 
to the uniformity of the distribution of the cpi. For a more complete characterization of the 
statistics of the traces we propose to calculate their marginal probability densities 




(29) 



where 5 2 (t) is a two dimensional delta function in the complex t plane; the ensemble average is 
to be done with the weight (0) . It turns out convenient to first calculate the Fourier transform 



PN.n( k ) = eX P 



-V 

2 ^ 



(30) 



Due to the periodicity of all functions of the phases ifi involved the characteristic function 



Ppi (k) depends on k only through the modulus |fc| while the density Pjy n {t) is only a function 



of \t\. Henceforth we assume k — k* — \k\ and write 



P Nn( k ) = ex P -ifc V"cos(n<^i) 



(31) 



We immediately obtain for the Poissonian ensemble 
P° N (k) = J (k) N , 



(32) 



Jo(fc) being a Bessel function. It follows that the densities P^ are the same for all values of 
the exponent n, a rather intuitive result given the statistical independence of the phases (pi in 
the Poissonian case. Equally expected for such a Poissonian random walk is the independence 
of the mean squared "displacement" of n, \t n \ 2 = N. 

For the unitary ensemble, (3 = 2, we exploit the symmetry of the exponential to be averaged 
in the phases (fii and employ the analogue of (^|) and find the characteristic function to take the 
form of a Toeplitz determinant, 



dip 
2^ 



The remaining <p integral again yields Bessel functions and can be written as 
l (Jb)=det| >^ J| s |(fc)(-i) |s| 5(£-m + ns) 



N. 



p2 
N. 



/ +oo 

det Y. J M( fc )(-i) 

\S — — OQ 



(33) 



(34) 



Clearly, the number of non-zero elements of the N x N determinant here incurred decreases as 
the order n of the trace t n in consideration grows. In particular, for n > N only the diagonal 
elements are non-zero such that 



Jo(k) 



N 



for n > N. 



(35) 
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At the other extreme, n = 1, we meet with the full Toeplitz determinant 



p2 



i(k) = T J 



N 



Jo 


-iJi 


-J 2 


-iJl 


Jo 


-iJi 


-Jl 


-iJi 


Jo 



(36) 



Intermediate values of n lead to the subdeterminants T m obtained by cancelling the last 
N — m rows and columns of TV, 



N-n 



■x, 



2n-N 



P^Jk)=T 3 N -^-T : 



n3n-N 



N 

for N > n> — , 
- - 2 ' 

N N 
tor — > n > — 
2 - - 3 



(37) 



and so forth. This can be seen as follows. Starting with N > n > ^ one checks that the 
determinant in ( |34| ) has non- vanishing elements residing only in the diagonal, Jo(k)5(£ — m), 
and in two subdiagonals, — iJi(k)5(£— m+n). One moves the (l+n)th row to become the second, 
then the (l+n)th column to become the second and thus isolates a 2x 2 block Ti = Jo(fc) 2 + Ji(fc) 
in the upper left corner. One repeats this process by moving the (2 + n)th rows and columns to 
become the fourth and so forth until one arrives at a block diagonal determinant in which the 
2x2 block T2 and the lxl block T\ = Jo(k) appear N — n times and 2n — N times, respectively. 
The procedure for ^ > n > y is analogous: One moves the (1 + n)th and the (2 + rt)th row 
to become the second and third, respectively, then does likewise to the (1 + n)th and (2 + n)th 
column and thus generates the 3x3 block T3 and so forth. As n decreases towards 1 we meet 
all the 



p2 (i\ _ rpN—mn rp(m+l)n- 
r N,n\ h -) — 1 m+l ' 1 m 



N 



N N 

for — > n > 

m m+1 



(38) 



with m = 1, . . . N. 

Simplest to deal with is, of course, the case of the smallest non-trivial dimension, N 
and there we obtain 



P?,i(k) = J%{k) + JfW, Pl 2 {k) = J 2 {k). 
By Fourier transforming we produce the densities of the first and second trace, 



\t\ 2 



27T 2 |t| 



P 2 



2:2 



(*) 



(39) 



(40) 



% n of the nth trace behaves like a Gaussian 



Now we propose to show that the distribution P l 

1 2m f or sufficiently low orders. To this end we define the functions 



with respect to its moments \t r 



J rn (k) = (-ik) m J m (k) 
which have the property (Abramowitz and Stegun 9.1.30 pl|) 
i d 



k dk 



Jm (k) = Jm-1 (k) . 



The Toeplitz determinant Tjy can then be rewritten as 

T N = det (J m , J m -i, .. .J , with m = 0, 
From this we find with (|lj) and @ for I < N 



,N-1. 



k=Q — 



-1 



(41) 



(42) 



(43) 



(44) 



To prove the foregoing identity we proceed as follows. Applying | ^ once to Tjv we get a 
determinant differing from Tjv only in the last column where according to ( |4^ ) J m -N+i—>Jm-N- 
Now setting k = and invoking J m (0) = for positive integer n while J m (0) = (^ i ) m /(— m)! 
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for m = 0, — 1, —2, . . . we face a determinant with vanishing elements below the diagonal and 

thus equalling the product of its diagonal elements; among these only the last, Ji (0) = — i/2, 
differs from unity whereupon ( f4~4| ) is proven for 1=1. In the next step, I = 2, TV is changed 
such that only the elements in the last 2 columns may get their indices shifted; again setting 



k = we are left with a 2 x 2 determinant which yields (44) for 1 — 2. In the ^th step we 
get a sum of I x I determinants which is evidently independent of the dimension N for I < N. 
Due to that independence we need not pursue the non-trivial task of proving ( fiij ) for arbitrary 
i < N < oo but rather invoke the much more easily proven result, to be established in the next 
section, that for N — > oo the determinant T/v is a Gaussian in A:. It follows that for finite ./V 
the expansion of T/v in powers of k 2 coincides with that of exp f — if) U P ^° Nth. power. 
Equations ([55]) - (37) show that n (t) behaves like a Gaussian distribution with respect to all 
moments |£„| 2m of orders m < N/n. In particular, we find for the variances of the traces 

IFF- J n fOT 0<n<iV 

IM -\ N for n>N. ^ 4&j 

Interestingly, these variances grow towards the Poissonian value N as n — > N from below and 
then remain stuck as n grows further. 

We now proceed to the orthogonal and symplectic cases. Starting as in the previous section 
with integrating over alternating variables we find for the orthogonal case with even N and the 
symplectic case 

Pk,„(k) oc (detA) 1 / 2 for 0=1 

PN.n{ k ) « (detB) 1/2 for (5 = 4 (46) 

with the antisymmetric matrices 

(-i)l s l(-i)l s 'l N-l 
A mm > = 2^ J\ s \{k)J\ s '\{k) — 8 (m + m + n(s + s )) with \m\, \m \ < — - — 

ss' 

B mml = J2 J \s\( 2 k)(-i) ls] 6(m + m' + ns) with \m\,\m'\ < 2 ~ 1 . (47) 

s 

We have not proven the Gaussian property but have calculated the variances. In the orthogonal 
case (j3 = 1) we find 

f 2n -n JT , tt~Jf ttto for < n < N 

2 AT — n > , — ; A. , -, Wn lor n> N , 

{ L-im=\ m+n— (JV+l)/2 — ' 

while the symplectic case (/? = 4) yields 

I Jv for n > 2N . 

Needless to say, these mean squared traces could have been read off the well-known two-level 
correlation functions of the circular ensembles the Fourier transforms of which our variances in 
essence are Indeed, by introducing a non- normalized density of eigenphases as 

N N oo 

P(V») = 2tt£;%- W )=X; E e~ in ^-^ (50) 

i—1 i—1 n=—oo 

one immediately sees that the two-point correlation function of that density reads 

oo 

P(<P)P(V')= E P^'- (51) 



n— — oo 



Figure ^ reveals excellent agreement of the mean squared traces as computed for samples of 
4 x 10 4 20 x 20 matrices of the four ensembles considered with the analytical predictions. Note 
that for the Poisson circular ensemble \t n \ 2 equals the matrix size N, without dependence on 
n. For small n the variance of traces decreases with the repulsion parameter /3. The data for 
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Figure 4. Mean squared traces \t n \ 2 for 10 4 matrices of size N = 20 typical 
(□), CUE (O ) and CSE (A) compared with analytical results (line! 
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Figure 5. Distribution of the modulus of the first trace P(|ti|) for 10 5 CUE matrices for a) 
N = 2 (semicircle distribution), and b) TV = 20 (Gaussian distribution). 
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Figure 6. Distribution of the moduli of the first five traces ti|, . . . , \ t$\ for 10 5 CUE matrices 
of size TV = 20. Narrower solid lines correspond to Gaussian distributions with appropriate 
variances \t n \ 2 = n. 



p2 



the symplectic ensemble are obtained with 2N dimensional matrices, which provide N different 
eigenvalues each. 

Figures [|, [| display the similarly fine agreement of our analytical results for the densities 
with numerical data for sets of random matrices drawn from the circular unitary ensemble. 
We have performed similar numerical studies for the orthogonal and symplectic ensembles, again 
finding Gaussian marginal distributions for the traces of sufficiently large matrices. Moreover, in 
all cases studied, the distribution of the traces was isotropic, i.e. without any phase dependence. 

We should add a word of intuitive explanation to the statistics of the traces for large 
dimensions N. The nth trace of a unitary matrix may be thought of as a random walk in 
the complex plane, with each of the N steps of unit length and the ith step in a direction 
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given as n times the ith eigenphase ipi. These directions are mutually independent and 
uniformly distributed for the Poissonian ensemble which fact explains the independence of the 
characteristic function ( |32] ) of n and the ensuing first and second moments, t n = 0, |t„| 2 = N . 
In accordance with the central limit theorem the rescaled traces In/^/N tend, with N — > oo, 
to have a Gaussian distribution of zero mean and unit variance. Without such rescaling the 
moments |i n | 2m with m <C N differ from those of the Gaussian defined by the first two moments 
only by corrections of relative order 1/N. 

For the circular ensembles with (3 > the phases ipi display repulsion of degree (3 such that 
the directions of subsequent steps in the random walk mentioned are not independent. The 
correlations between the phases cannot prevent near-Gaussian behaviour of the traces t n with 
n <C N, as is intuitive in view of the local character of the spectral correlations. Moreover, 
while the phases tpi cover the interval [0, 2n) uniformly once, their multiples rupi go around 
that interval n times such that for n > 1 the phases nifi [mod (2ir)] may exhibit accidental close 
neighborhoods of originally distant <pi. 



4. Joint density of traces for large CUE matrices 

We shall here employ a powerful theorem about determinants of Tocplitz matrices, due originally 
to Szego and Kac and extended by Hartwig and Fisher p2| , to find the marginal and joint 
distributions of the traces t n of CUE matrices in the limit of large dimension N. As our starting 
point we recall the identity (||) for the CUE average of a symmetric function of all N phases. 
Assuming, moreover, that symmetric function to have the form of a product we can pull the 
integral over the mth phase tp m into the mth row of the determinant in (Q) and thus express 
the average as a Toeplitz determinant, 



N 



Y[ f{<p m ) = det(/ m , f m _ u f- N+1 ) = T ({/}) , m = 0,l,...,7V-l, (52) 



m—l 

the elements of which are given by the Fourier transform 

r 2ir 



(53) 



of the function f{f). We had incurred two examples in ([|), ( |io| ) and (|33|), (|3^). The theorem in 
question says that for large N the above determinant is given by 

oo 

In T ({/}) = Nl + nl n l- n (54) 

n=0 

where the l n are the Fourier coefficients of In /(</?), i.e. In f(ip) = X)^L-oo ^n^ mip ■ The conditions 
the function f(<p) must meet for the above limiting form to hold are (i) f(ip) for < <p < 2tt, 
(ii) arg/(27r) = arg/(0), (iii) E,T=-oo \fn\ < and ( iv ) E^-oo l"-ll/"l 2 < °°; the y are fulfilled 
in all examples of interest here. 

In a first application we return to the Toeplitz determinant ( [34] ) which gives the density 
Pjq n of the nth trace. The only non- vanishing Fourier coefficients of ln/(<^) = — ikcosrup are 
l± n = —ifc/2 whereupon we get the Gaussian anticipated in the previous section, 

P£ ;n (fc) = e-« fe2 / 4 ^ Pl n (t) = — e-l*! 2 /", for N » 1. (55) 

rm 

No more difficult to obtain is the joint density of the first n traces 



m=l \ i=l / 

since its Fourier transform P^(ki, . . . , k n ) is once more of the form (^2|) with 

f(<p) = ex P ( -i ( k ^ imV + K n e imv )) ■ (57) 

\ m=l / 
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The non-vanishing Fourier coefficients of the logarithm of that latter function are l r , 
—ik m /2, l- m = — \k* m j2 with m — 1, . . . , n. The theorem (|54| ) thus yields, for N 3> n, 

P^(ki,...,k n ) = exp (- m\k m \ 2 /4\ 



^Pff(t u ...,t n ) = _expl - ^ MVmj, (58) 

\ m=l / 

i. e. the product of the marginal distributions of the first n traces. The result generalizes in an 
obvious way to the joint density of an arbitrary set of finite-order traces. We thus conclude that 
in the limit N — > oo the finite-order traces are statistically independent and all have Gaussian 
distributions. 

We can now briefly comment on the conditions of applicability of the Hartwig-Fisher theorem 
given above. The first two of them are clearly fulfilled here since i In /(<£>) as given by (EtI) 
is real, continuous and periodic. The third and fourth condition are met since the derivative 
f'if) — Em=-oo lm fmZ lm ' p is square integrable (trivially indeed since iln/(y>) is a finite Fourier 
series); in particular, f'(ip) obeys Parseval's identity, 

oo 2 5 n 

E M 2 |/J 2 = / |/'M| 2 d^ = 7rEm 2 |fc m | 2 <oo. (59) 

m— — oo q m— 1 

Since |m||/ m | 2 < |m| 2 |/ m | 2 it follows that X!m=-oo l m ll/m| 2 < oo, i-e* the validity of condition 
(iii)- 

Finally, we invoke Cauchy's inequality in 



L l/™l = l/ol+L 

m= — oo m^O 



\mf m \ < |/o 





(60) 



Upon using (|59|) and the convergence of X^m^o p^Tp we ver ify condition (iv). 

For finite dimension the independence as well as the Gaussian character of the traces are 
only approximate. For sets of traces both of these properties tend to get lost as the sum of the 
orders of the traces in a set increases. In particular, since all traces t n are uniquely determined 
by the N real eigenphases, only N/2 traces can be independent. 

Preliminary numerical studies suggest that the finite-order traces might be similarly 
independent and Gaussian for the COE and the CSE. For the CPE, of course, the independence 
holds trivially. 



5. Remarks on the distributions of coefficients of secular polynomials 

As we have just seen the first few traces 1 1 , ti , . . . of large unitary matrices drawn from any of the 
circular ensembles display no noticeable correlations. It follows that the coefficients a\ , a-i , . . . 
must bear strong mutual correlations simply since a n can be expressed in terms of the first n 



traces through Newton's formulae (25). One must therefore expect that the explicit form of the 
marginal and joint distributions of the a n are hard to come by. An exception is provided by a\ 
the marginal distribution of which is trivially related to that of the first trace since a\ =t\. 

A slightly less trivial result may be obtained for the distribution P(a,2) of the second 
coefficient: Since ai = \(t\ — t-i) we may invoke the CUE joint distribution (58) of the first two 
traces to get 

Pfa) = a/5 f r^-^ v (-- ^ 2x ) ■ ^ 

V 7T 3 J Q y/ X (2-x) V X J 

A saddle-point approximation to the foregoing integral immediately reveals that P{a,2) decays 
exponentially for large a 2 . Proceeding similarly one may combine Newton's formulae with the 
joint distribution of the traces to get the marginal and joint distributions of the first few a„, 
with decreasingly compact and enjoyable results. 

Figure presents the distribution of the modulus of |ot2 1 as obtained numerically from 10 5 
CUE matrices of the size 10 x 10. These numerical data agree well with the distribution 
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1.0 




Figure 7. Distribution of the moduli of second coefficients 1 0.2 1 for 10 CUE matrices of size 
N = 20. Solid line represents 2w\a2\P(a2) according to (5.1) and is compared to a dashed line 
representing Gaussian distribution of |ai|. 




Figure 8. Distribution of the moduli of the coefficients |ai| (short dashed line), |<X2 1 (dotted), 
\<14,\ (dash-dotted) and \a&\ (long dashed line) for CUE matrices of size N = 20. Distributions of 
the two first coefficients are compared with theoretical predictions (narrow lines). 



2-!r\a2\P(ci2) according to (|6l|). In Figure || we depict numerical results for the distributions 
of the moduli of some higher coefficients \a n \. All of these curves grow linearly out of the origin 
with a slope increasing with the index n of the coefficient a n . On the other hand, all of these 
distributions but the first are characterized by a long exponential tail. The latter originates 
from the convolution type integrals which combine the densities of the traces to those of the 
coefficients; it contrasts with the Gaussian tail of the distribution of the first coefficient. The 
qualitative features just outlined for the CUE are shared by the distributions P(a n ) for the 
other canonical ensembles. 



6. Comparison with a dynamical system 

We here propose to examine in how far our results on distributions of the coefficients a n and 
traces t n of matrices drawn from the canonical circular ensembles are applicable to dynamical 
systems. Choosing the familiar model of the periodically kicked top ||, (23) we work with a 
Floquet operator of the form 

F = cxp (y-' l 2j X + \ J * ~ ipx exp (~ i 2j~+T J ^ ~ [Pv Jy 

x exp f -i-A- - i Pz J z ) . (62) 



2j + l 

This involves the components of an angular momentum operator J Xl J yi Jz which satisfy 
standard commutation relations, [J x , J y ] = iJ z etc. The quantum number j fixes J 2 = j(j + 1) 
and the size of the Hilbert space, N = 2j + 1. 

For generic values of the parameters k x ,k y ,k z and p x ,p y ,p z the corresponding classical 
dynamics is chaotic and there is no geometric nor antiunitary symmetry left |^| . All previously 
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studied statistical properties of the quasienergy spectrum and the eigenvectors were found to be 
remarkably faithful to the predictions of the CUE H . On the other hand, when the parameters 
k y and p y (or, instead, k x and p x ) are set to zero, an antiunitary symmetry under time reversal 
appears, and in that case the spectral and eigenvector statistics were found as of the COE type. 
We shall refer to the two variants of the model as to the unitary and orthogonal top, respectively. 

Before presenting our data for the traces t n and the secular coefficients a n of various tops a 
word of caution is in order. Previous statistical analyses of tops were made for spectrally local 
quantities like low-order correlations of the level density, the distribution of nearest-neighbor 
spacings, or for components of eigenvectors; what distinguishes these quantities is a certain self- 
averaging character: A single Floquet matrix of large dimension N provides a sufficient data 
basis to extract reliable means or even distributions from. Now, the sequence of the traces t n 
and that of the secular coefficients a n are not in any way self-averaging since such a sequence 
with 1 < n < N/2 stands for a whole quasienergy spectrum and therefore changes in a system 
specific manner when control parameters are varied. Consequently, it would not make sense to 
compare such a sequence for an individual Floquet matrix with the means calculated here for 
the various circular ensembles. We must rather ask whether an ensemble of Floquet matrices of 
the type (62), defined by a whole set of values for the various control parameters, is faithful to 
the prediction based on the circular ensemble of random matrices of the same symmetry class. 
It is in this sense that we have undertaken the comparison to follow. 




n 



Figure 9. Mean squared traces \t n \ 2 (circles) and coefficients \a n \ 2 (triangles) for orthogonal 
tops are compared with COE results (lines). Data are averaged over 20000 different Floquet 
matrices of the structure (6.1) as described in the text. 

Figure || shows the n-dependence of the variances of the traces (circles) and of the coefficients 
(triangles) for orthogonal tops. The order of the characteristic polynomial was taken as N = 41 
by choosing j = 20. Data from 20000 matrices were gathered by picking k y , k z from intervals of 
length 3 around 10 and p yi p z from [3tt/8, 57r/8], all with independent box distributions. These 
intervals were chosen so as to secure classical chaos and to avoid geometric symmetries. The 
agreement with the COE is obviously satisfactory. 
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kicked top from CSE 




5 10 15 20 25 30 35 40 



n 

Figure 11. As in Figure ^ for symplectic tops. Average over 2000 matrices as described in the 
text. Lines denote CSE results. 



Analogous data for unitary tops are presented in Figure |l0|. The agreement with the 
CUE predictions is acceptable in the sense that there is no doubt about the universality class. 
However, system specific behaviour is clearly visible for low-order traces and even in all secular 
coefficients. Qualitatively, such deviations from random-matrix theory are not unexpected since 
low-order traces do contain primarily system specific information retrievable from short periodic 
orbits and since the a n even for large n contain low-order traces as expressed in Newton's 
formulae. 

Symplectic tops can also be constructed by securing an antiunitary symmetry T with 
T 2 = — 1 and avoiding geometric symmetries jjj . The Floquet operator 

F = exp ( _i ^TT Jl i WTT) {Jx Jz + JzJx) ~ i WTT) {Jx Jy + JyJx) ) 

x exp (-i^f) (63) 

in a representation with halfinteger j serves fine. A set of 2000 matrices with j = 39.5 and 
coupling constants k{ drawn at random from a hypercube of length 0.4 near 10 gives data 
shown in Figure [ll]. The agreement with the CSE is better than could be hoped for. 
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Appendix A. Differential equation for the generating function 



Let us recall the heuristic way we first went towards the mean squared secular coefficients (19), 
(p3|), and (24). After deriving the Pfaffians ( |l8| ) and ( pi] ) we applied an algebraic computer 
program to evaluate the \a n \ 2 for small dimensions N and then proceeded analytically to n — 1,2 
and arbitrary N. We thank N Lehmann for the computer work and also for his intuition in 
guessing ( p3|) with us. Having eventually guessed ( |l9| ) as well we speculatively extrapolated to 
the general form (24) for arbitrary positive (3. 

To prepare for the proof of (24) we observe the equivalence to the recursion relation 

(n + 1) (l + —(N — n— I)) = (N-n)(l + n^j (Al) 



with | a,o | 2 = 0. Multiplying (Al) by x n , using (x-f^)x n = nx n , and summing from n = 
to n = N — 1 we obtain the differential equation (|25|) the polynomial solution of which with 
Pj^(x) = 1 generates the mean squared secular coefficients. A few elementary steps will now 
yield the validity of the differential equation. 
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We start from the definition (4) of the generating function, setting A = x, /i — 1 and take 
the derivative with respect to x, 



dP(x) 
dx 



N 1 



(A2) 



where the angular brackets denote an average with the whole integrand of ( |A2|) as the weight. 
On the other hand, taking the derivative with respect to /i before setting /i = 1 and exploiting 
the invariance of the A-fold phase integral under a constant shift of all phases we get 



dP(x) 
dx 



^ e^* - 1 



(A3) 



Next, we differentiate (A2) once more 
d 2 P{x) 



dx 2 



y I 

(pi 



(e^> - x)(e^i - ew) 



y- 

Z-^i (pi 



-llfii 



(e 1 ^ - x) 



1 — i cot 



<Pi - V] 



(A4) 



Realizing that the cotangent function may be replaced by a derivative acting on the joint density 
of eigenvalues (0) , 



9y> 



(A5) 



we can transform the last member in (A4) by partial integration. A little algebraic hocuspocus 
then gives, with the help of (A2) and (A3), 



2 P(x) 
dx 2 



A- 1 + 



0{ 



^-x) \l_(dPy\ 

0(1 - x) ) x yox y / 



(A6) 



Finally, the average J2i( e can be expressed in terms of the generating function P and its 
derivative dP/dx by employing the identity X]fe( e_IV,fc ) = J2k(^'^: e ~ liPk ) > integrating by parts, 
and again invoking ( A2 , A|, A5 ) , 



£ < e -Vi) ( 2+ |(JV -!))=(*-!) 



dP 

dx 



NP. 



(A7) 



Upon inserting this in (A6) we obtain a differential equation for the generating function P which 
is easily put into the form (25). 
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